7-3 all cause mortality

Analyses of all cause mortality.

Authors
Affiliations

James Totterdell

University of Sydney

Rob Mahar

University of Melbourne

Published

July 26, 2023

Load packages
library(ASCOTr)
library(tidyverse)
library(lubridate)
library(kableExtra)
library(patchwork)
library(cmdstanr)
library(posterior)
library(bayesplot)
library(ggdist)
library(lme4)
library(broom)
library(broom.mixed)
library(bayestestR)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))

bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))

color_scheme_set("red")
options(digits = 4)
Prepare datasets
all_dat <- readRDS(file.path(ASCOT_DATA, "all_data.rds"))

# FAS-ITT
fas_itt_dat <- ASCOTr:::make_fas_itt_set(all_dat)
fas_itt_nona_dat <- fas_itt_dat |>
  filter(!is.na(D28_death))

# ACS-ITT
acs_itt_dat <- ASCOTr:::make_acs_itt_set(all_dat)
acs_itt_nona_dat <- acs_itt_dat |>
  filter(!is.na(D28_death))

# AVS-ITT
avs_itt_dat <- ASCOTr:::make_avs_itt_set(all_dat)
avs_itt_nona_dat <- avs_itt_dat |>
  filter(!is.na(D28_death))
Load models
logistic_mod <- compile_cmdstanr_mod(
  file.path("binary", "logistic"), dir = "stan")
logisitc_site <- compile_cmdstanr_mod(
  file.path("binary", "logistic_site"), dir = "stan")
logistic_site_epoch <- compile_cmdstanr_mod(
  file.path("binary", "logistic_site_epoch"), dir = "stan")
Code
make_D28_death_table_A <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(AAssignment = factor(
      AAssignment, 
      levels = paste0("A", 0:2),
      labels = intervention_labels()$AAssignment)
    ) %>%
    group_by(AAssignment) %>%
    summarise(
      Randomised = n(),
      `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(D28_death)), 100 * mean(is.na(D28_death))),
      `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(D28_death)), 100 * mean(!is.na(D28_death))),
      `Died within 28 days` = sprintf("%i (%.1f)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(AAssignment = "Overall") %>%
        summarise(
          Randomised = n(),
          `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(D28_death)), 100 * mean(is.na(D28_death))),
          `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(D28_death)), 100 * mean(!is.na(D28_death))),
          `Died within 28 days` = sprintf("%i (%.1f)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
      )
    ) %>%
    mutate(AAssignment = fct_inorder(AAssignment)) %>%
    gather(key, value, -AAssignment, factor_key = T) %>%
    spread(AAssignment, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(
      font_size = 9, 
      latex_options = "HOLD_position")  
}


make_D28_death_table_C <- function(dat, format = "html") {
  tab <- dat %>%
    mutate(CAssignment = factor(
      CAssignment, 
      levels = paste0("C", 0:4),
      labels = intervention_labels()$CAssignment)
    ) %>%
    group_by(CAssignment) %>%
    summarise(
      Randomised = n(),
      `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(D28_death)), 100 * mean(is.na(D28_death))),
      `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(D28_death)), 100 * mean(!is.na(D28_death))),
      `Died within 28 days` = sprintf("%i (%.1f)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
    ) %>%
    bind_rows(
      dat  %>%
        group_by(CAssignment = "Overall") %>%
        summarise(
          Randomised = n(),
          `Outcome missing` = sprintf("%i (%.1f)", sum(is.na(D28_death)), 100 * mean(is.na(D28_death))),
          `Outcome observed` = sprintf("%i (%.1f)", sum(!is.na(D28_death)), 100 * mean(!is.na(D28_death))),
          `Died within 28 days` = sprintf("%i (%.1f)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
      )
    ) %>%
    mutate(CAssignment = fct_inorder(CAssignment)) %>%
    gather(key, value, -CAssignment, factor_key = T) %>%
    spread(CAssignment, value)
  colnames(tab)[1] <- "n (\\%)"
  if(format == "latex") {
    colnames(tab) <- linebreak(colnames(tab), align = "c", linebreaker = "<br>")
  }
    kable(tab,
      format = format,
      align = "lrrrrrr",
      escape = FALSE,
      booktabs = TRUE
    ) %>%
    kable_styling(
      font_size = 9, 
      latex_options = "HOLD_position")  
}

Descriptive

Anticoagulation

Relationship anti-coagulation to outcome
make_D28_death_table_C(
    acs_itt_dat
)
n (\%) Low
dose
Intermediate
dose
Low dose
with aspirin
Therapeutic
dose
Overall
Randomised 610 613 283 50 1556
Outcome missing 12 (2.0) 5 (0.8) 2 (0.7) 0 (0.0) 19 (1.2)
Outcome observed 598 (98.0) 608 (99.2) 281 (99.3) 50 (100.0) 1537 (98.8)
Died within 28 days 19 (3.2) 15 (2.5) 10 (3.6) 6 (12.0) 50 (3.3)
Table 1: Primary outcome by anti-coagulation intervention.
Relationship anti-coagulation to outcome
save_tex_table(make_D28_death_table_C(
      acs_itt_dat,
      "latex"), 
      "outcomes/secondary/7-3-anticoagulation-summary")

Antiviral

Relationship anti-viral to outcome
make_D28_death_table_A(
     avs_itt_dat
)
n (\%) Standard
of care
Nafamostat
Overall
Randomised 73 82 155
Outcome missing 0 (0.0) 0 (0.0) 0 (0.0)
Outcome observed 73 (100.0) 82 (100.0) 155 (100.0)
Died within 28 days 0 (0.0) 0 (0.0) 0 (0.0)
Table 2: Primary outcome by anti-viral intervention.
Relationship anti-viral to outcome
save_tex_table(make_D28_death_table_A(
     avs_itt_dat,
      "latex"), 
      "outcomes/secondary/7-3-antiviral-summary")

Age

Relationship age to outcome
agedat <- fas_itt_dat %>%
  dplyr::count(D28_death, AgeAtEntry) %>% 
  spread(D28_death, n, fill = 0) %>% 
  mutate(
    n = `0` + `1` + `<NA>`,
    p = `1` / (`1` + `0`))
agemod <- glm(
  cbind(`1`, `0`) ~ AgeAtEntry, 
  data = agedat, 
  family = binomial())
agedat <- agedat %>%
  mutate(
    ypred = predict(agemod, newdata = agedat, type = "response")
  )
p1 <- ggplot(agedat, aes(AgeAtEntry, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry")
p2 <- ggplot(agedat, aes(AgeAtEntry, p)) +
    geom_point() +
    geom_vline(xintercept = 60, linetype = 2) +
    geom_line(aes(y = ypred)) +
    labs(y = "Proportion\ndied by day 28", x = "Age at entry")
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-age.pdf"), p, height = 2, width = 6)
p

Figure 1: Relationship (logistic regression linear in age) between age at entry and the primary outcome.

Sex

Relationship sex to outcome
tdat <- all_dat %>%
  filter_fas_itt() %>%
  dplyr::count(Sex, D28_death) %>%
  group_by(Sex) %>%
  spread(D28_death, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(Sex, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants")
p2 <- ggplot(tdat, aes(Sex, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\ndied by day 28")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-sex.pdf"), p, height = 2, width = 6)
p

Figure 2: Mortality by sex.

Oxygen

Relationship oxygen to outcome
tdat <- all_dat %>%
  filter_fas_itt() %>%
  dplyr::count(
    supp_oxy2 = factor(supp_oxy2, labels = c("Not required", "Required")), 
    D28_death) %>%
  group_by(supp_oxy2) %>%
  spread(D28_death, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(supp_oxy2, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", x = "Supplemental oxygen")
p2 <- ggplot(tdat, aes(supp_oxy2, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\ndied by day 28", x = "Supplemental oxygen")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-oxygen.pdf"), p, height = 2, width = 6)
p

Figure 3: Mortality by oxygen.

Country

Relationship country to outcome
tdat <- all_dat %>%
  filter_fas_itt() %>%
  dplyr::count(Country = fct_infreq(PT_CountryName), D28_death) %>%
  group_by(Country) %>%
  spread(D28_death, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  )
p1 <- ggplot(tdat, aes(Country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")
p2 <- ggplot(tdat, aes(Country, p_1)) +
  geom_point() +
    labs(
      y = "Proportion\ndied by day 28", 
      x = "Country of enrolment")  +
  ylim(0, NA)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-country.pdf"), p, height = 2, width = 6)
p

Figure 4: Mortality by country.

Site

Relationship site to outcome
tdat <- all_dat %>%
  filter_fas_itt() %>%
  dplyr::count(
    Country_lab = Country,
    Site_lab = fct_infreq(Location),
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand")),
    Site = PT_LocationName,
    D28_death) %>%
  group_by(Country, Site) %>%
  spread(D28_death, n, fill = 0) %>%
  mutate(
    n = `1` + `0` + `<NA>`,
    p_1 = `1` / (`1` + `0`),
    p_miss = `<NA>` / (`1` + `0` + `<NA>`)
  ) %>%
  ungroup()
p1 <- ggplot(tdat, aes(Site_lab, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(tdat, aes(Site_lab, p_1)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_point() +
    labs(
      y = "Proportion\ndied by day 28", 
      x = "Site of enrolment") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p <- p1 / p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-country-site.pdf"), p, height = 4, width = 6.25)
p

Figure 5: Day 28 mortality by site within country.

Calendar Time

Relationship calendar date to outcome
caldat <- all_dat %>% 
  filter_fas_itt() %>%
  dplyr::count(D28_death, yr = year(RandDate), mth = month(RandDate)) %>% 
  spread(D28_death, n, fill = 0) %>% 
  mutate(p = `1` / (`1` + `0`),
         n = `1` + `0` + `<NA>`)
p1 <- ggplot(caldat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_point() +
    labs(
      y = "Proportion\ndied by day 28", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12) +
  ylim(0, 0.15)
p2 <- ggplot(caldat, aes(mth, n)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p <- p2 | p1
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-3-calendar-time.pdf"), p, height = 2, width = 6)
p

Figure 6: Relationship between calendar date and the primary outcome.

Analyses

FAS-ITT

Code
res <- fit_primary_model(fas_itt_nona_dat, logistic_site_epoch, ctr = contr.equalprior, outcome = "D28_death")
res$drws$OR <- res$drws$OR[-1]
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Ineligible aspirin", "Age \u2265 60", "Female", "Oxygen requirement", "Australia/New Zealand", "Nepal")
Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-3-primary-model-fas-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Nafamostat 0.88 (0.10, 8.22) 1.68 (2.79) 0.55
Intermediate-dose 0.79 (0.39, 1.60) 0.85 (0.32) 0.74
Low-dose with aspirin 0.80 (0.34, 1.82) 0.87 (0.39) 0.70
Therapeutic-dose 3.95 (1.14, 13.02) 4.75 (3.26) 0.02
Ineligible aspirin 4.69 (1.15, 16.80) 5.79 (4.24) 0.02
Age ≥ 60 2.00 (1.05, 3.78) 2.10 (0.70) 0.02
Female 0.34 (0.16, 0.69) 0.37 (0.14) 1.00
Oxygen requirement 3.23 (1.67, 6.25) 3.42 (1.19) 0.00
Australia/New Zealand 0.53 (0.11, 2.36) 0.71 (0.62) 0.80
Nepal 2.51 (0.55, 9.80) 3.17 (2.52) 0.11
Posterior summaries for model parameters (fixed-effects), FAS-ITT.
Odds ratio summary for epoch and site
p <- plot_epoch_site_terms(
  res$drws$gamma_epoch,
  res$drws$gamma_site,
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-3-primary-model-epoch-site-terms-fas-itt.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Summary of epoch and site posterior odds ratios.
Code
p <- plot_or_densities(c(res$drws$AOR, res$drws$COR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-3-primary-model-fas-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Posterior densities for treatment comparisons.
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8726 0.8759 0.8284 0.8378 0.9078 0.8129 0.8534 0.8569
Code
mcmc_trace(res$drws["beta"])

Code
mcmc_trace(res$drws["gamma_site"])

Code
mcmc_trace(res$drws["gamma_epoch"])

Posterior predictive

Code
y_ppc <- res$drws$y_ppc
ppc_dat <- bind_cols(fas_itt_nona_dat, tibble(y_ppc = y_ppc))
grp_ppc <- function(grp) {
  ppc_dat %>%
    group_by(grp = !!grp) %>%
    summarise(
      mean_y = mean(D28_death),
      rvar_mean_y_ppc = rvar_mean(y_ppc)
    )
}
plot_grp_ppc <- function(dat, lab = "") {
  dat %>%
    ggplot(aes(y = grp, xdist = rvar_mean_y_ppc)) +
    stat_interval(size = 2) +
    geom_point(aes(x = mean_y, y = 1:nrow(dat)), colour = "red", shape = 23) +
    labs(
      x = "Posterior predictive\nproportion", 
      y = lab)  
}
ppc_A <- grp_ppc(sym("AAssignment"))
ppc_C <- grp_ppc(sym("CAssignment"))
ppc_ctry <- grp_ppc(sym("Country"))
ppc_epoch <- grp_ppc(sym("epoch"))
ppc_site <- ppc_dat %>%
  group_by(Country = Country, grp = site) %>%
  summarise(
    mean_y = mean(D28_death),
    rvar_mean_y_ppc = rvar_mean(y_ppc)
  )
p0 <- plot_grp_ppc(ppc_A, "Antiviral\nintervention") + labs(x = "")
p1 <- plot_grp_ppc(ppc_C, "Anticoagulation\nintervention") + labs(x = "")
p2 <- plot_grp_ppc(ppc_ctry, "Country") + labs(x = "")
p3 <- plot_grp_ppc(ppc_epoch, "Epoch") + labs(x = "")
p4 <- plot_grp_ppc(ppc_site %>% filter(Country == "IN"), "Sites\nIndia")  + labs(x = "")
p5 <- plot_grp_ppc(ppc_site %>% filter(Country == "AU"), "Sites\nAustralia") + labs(x = "")
p6 <- plot_grp_ppc(ppc_site %>% filter(Country == "NP"), "Sites\nNepal")
p7 <- plot_grp_ppc(ppc_site %>% filter(Country == "NZ"), "Sites\nNZ")
p <- ((p3 | p0 / p1 / p2) / 
        ( (p4 | p5) / (p6 | p7) + 
            plot_layout(heights = c(3, 1)) ) ) +
  plot_layout(heights = c(1, 1.5), guides = "collect")
pth <- file.path("outputs", "figures", "outcomes", "secondary",
                 "7-3-primary-model-fas-itt-ppc.pdf")
ggsave(pth, p, width = 6, height = 5.5)
p

Posterior predictive distribution versus observed proportion (red diamond).

ACS-ITT

Fit primary model
res <- fit_primary_model(acs_itt_nona_dat, logistic_site_epoch, ctr = contr.equalprior, outcome = "D28_death")
res$drws$OR <- res$drws$OR[-1]
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Ineligible aspirin", "Age \u2265 60", "Female", "Oxygen requirement", "Australia/New Zealand", "Nepal")
Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-3-primary-model-acs-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Nafamostat 0.83 (0.08, 8.54) 1.68 (3.11) 0.56
Intermediate-dose 0.80 (0.39, 1.64) 0.85 (0.33) 0.73
Low-dose with aspirin 0.82 (0.35, 1.85) 0.89 (0.39) 0.68
Therapeutic-dose 4.00 (1.18, 13.20) 4.83 (3.31) 0.01
Ineligible aspirin 4.48 (1.15, 16.34) 5.57 (4.16) 0.02
Age ≥ 60 2.04 (1.08, 3.85) 2.15 (0.72) 0.01
Female 0.34 (0.16, 0.69) 0.36 (0.14) 1.00
Oxygen requirement 3.24 (1.68, 6.36) 3.44 (1.20) 0.00
Australia/New Zealand 0.54 (0.11, 2.44) 0.72 (0.64) 0.79
Nepal 2.55 (0.54, 10.10) 3.24 (2.65) 0.11
Posterior summaries for model parameters (fixed-effects), ACS-ITT.
Odds ratio summary for epoch and site
p <- plot_epoch_site_terms(
  res$drws$gamma_epoch,
  res$drws$gamma_site,
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-3-primary-model-epoch-site-terms-acs-itt.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Summary of epoch and site posterior odds ratios.
Code
p <- plot_or_densities(c(res$drws$AOR, res$drws$COR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-3-primary-model-acs-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Posterior densities for treatment comparisons.
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.8030 0.7854 0.8184 0.8520 0.8222 0.8149 0.7771 0.7611
Code
mcmc_trace(res$drws["beta"])

Code
mcmc_trace(res$drws["gamma_site"])

Code
mcmc_trace(res$drws["gamma_epoch"])

AVS-ITT

Pre-specified

  • excludes epoch
Code
res <- fit_primary_model(
  avs_itt_nona_dat, 
  logisitc_site, 
  ctr = contr.equalprior, 
  outcome = "D28_death", 
  vars = c("agegte60", "sexF", "supp_oxy2", "crp_tertile", "ctry"), 
  beta_sd_var = c(2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 1)
)
res$drws$OR <- res$drws$OR[-1]
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Age \u2265 60", "Female", "Oxygen requirement", "CRP (2nd tertile)", "CRP (3rd tertile)", "CRP (unknown)", "Nepal")
Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-3-primary-model-avs-itt-summary-table-pre-spec")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Nafamostat 0.96 (0.08, 11.25) 2.15 (4.91) 0.51
Intermediate-dose 1.00 (0.08, 12.46) 2.29 (4.91) 0.50
Low-dose with aspirin 1.22 (0.08, 16.99) 2.99 (6.36) 0.44
Therapeutic-dose 1.15 (0.08, 15.97) 2.85 (6.71) 0.46
Age ≥ 60 0.37 (0.00, 14.12) 2.06 (8.65) 0.69
Female 0.29 (0.00, 9.14) 1.37 (6.11) 0.75
Oxygen requirement 0.15 (0.00, 4.32) 0.64 (2.36) 0.86
CRP (2nd tertile) 0.34 (0.00, 12.50) 1.80 (7.03) 0.71
CRP (3rd tertile) 0.34 (0.00, 12.54) 1.92 (8.69) 0.71
CRP (unknown) 0.49 (0.01, 23.02) 3.38 (19.52) 0.63
Nepal 0.96 (0.14, 6.73) 1.56 (1.89) 0.52
Posterior summaries for model parameters (fixed-effects), AVS-ITT.
Code
p <- plot_or_densities(c(res$drws$AOR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-3-primary-model-avs-itt-odds-ratio-densities-pre-spec.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Posterior densities for treatment comparisons.
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9217 0.9981 0.9607 1.0189 0.9518 1.0084 1.0286 1.0321
Code
mcmc_trace(res$drws["beta"])

Code
mcmc_trace(res$drws["gamma_site"])

Code
mcmc_trace(res$drws["tau_site"])

Reduced Model

  • excludes epoch, region, and site
Code
res <- fit_primary_model(
  avs_itt_nona_dat, 
  logistic_mod, 
  ctr = contr.equalprior, 
  outcome = "D28_death", 
  vars = c("agegte60", "sexF", "supp_oxy2", "crp_tertile"), 
  beta_sd_var = c(2.5, 2.5, 2.5, 2.5, 2.5, 2.5)
)
res$drws$OR <- res$drws$OR[-1]
names(res$drws$AOR) <- "Nafamostat"
names(res$drws$COR) <- intervention_labels2()$CAssignment[-(1:2)]
names(res$drws$OR) <- c("Age \u2265 60", "Female", "Oxygen requirement", "CRP (2nd tertile)", "CRP (3rd tertile)", "CRP (unknown)")
Code
XX <- res$dat$X
XX_m <- colMeans(XX)
XX <- cbind(XX[, 1:6], sweep(XX[, 7:12], 2, XX_m[7:12]))
mdat <- res$dat
mdat$X <- XX
snk <- capture.output(
  mfit <- logistic_mod[["sample"]](
    data = mdat,
    seed = 1234,
    adapt_delta = 0.95,
    refresh = 0,
    iter_warmup = 1000,
    iter_sampling = 2500,
    chains = 4,
    parallel_chains = min(4, parallel::detectCores())
  )
)
mpars <- mfit$metadata()$model_params
keep <- mpars[!grepl("(_raw|epsilon_)", mpars)]
mdrws <- as_draws_rvars(mfit$draws(keep))
names(mdrws$beta) <- colnames(mdat$X)
mdrws$OR <- exp(mdrws$beta[!grepl("(Intercept|rand)", names(mdrws$beta))])
median(mdrws$beta)
Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR), "latex"),
  "outcomes/secondary/7-3-primary-model-avs-itt-summary-table")
odds_ratio_summary_table(c(res$drws$AOR, res$drws$COR, res$drws$OR))
Parameter Median 95% CrI Mean (SD) Pr(OR < 1)
Nafamostat 0.99 (0.09, 11.38) 2.11 (3.82) 0.50
Intermediate-dose 1.00 (0.08, 12.95) 2.33 (4.53) 0.50
Low-dose with aspirin 1.21 (0.09, 16.23) 2.88 (6.31) 0.45
Therapeutic-dose 1.16 (0.09, 15.85) 2.83 (5.98) 0.46
Age ≥ 60 0.38 (0.00, 14.50) 2.13 (8.86) 0.68
Female 0.30 (0.00, 8.97) 1.47 (9.25) 0.74
Oxygen requirement 0.16 (0.00, 4.59) 0.71 (4.17) 0.85
CRP (2nd tertile) 0.35 (0.00, 13.59) 1.95 (7.31) 0.71
CRP (3rd tertile) 0.34 (0.00, 12.07) 1.79 (7.00) 0.71
CRP (unknown) 0.50 (0.01, 22.88) 3.26 (16.63) 0.63
Posterior summaries for model parameters (fixed-effects), AVS-ITT.
Code
p <- plot_or_densities(c(res$drws$AOR)) +
  labs(x = "Odds ratio (log scale)", y = "Comparison")
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-3-primary-model-avs-itt-odds-ratio-densities.pdf")
ggsave(pth, p, width = 6, height = 2.5)
p

Posterior densities for treatment comparisons.
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.0087 0.9952 1.0113 1.0508 0.9884 1.0704 1.0109 1.0048
Code
mcmc_trace(res$drws["beta"])

End of script.